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ABSTRACT 

A fundamentally multidimensional convection scheme is described based on vector transient 
interpolation modelling rewritten in conservative control-volume form. Vector third-order 
upwinding is used as the basis of the algorithm; this automatically introduces important cross- 
difference terms that are absent from schemes using component- wise one-dimensional formu- 
las. Third-order phase accuracy is good; this is important for coarse-grid large-eddy or full 
simulation. Potential overshoots or undershoots are avoided by using a recently developed 
universal limiter. Higher order accuracy is obtained locally, where needed, by the cost-effective 
strategy of adaptive stencil expansion in a direction normal to each control-volume face; this is 
controlled by monitoring the absolute normal gradient and curvature across the face. Higher 
(than third) order cross- terms do not appear to be needed. Since the wider stencil is used only in 
isolated narrow regions (near discontinuities), extremely high (in this case, seventh) order 
accuracy can be achieved for little more than the cost of a globally third-order scheme. 

INTRODUCTION 

The authors have recently developed (reference 1) a cost-effective strategy for obtaining 
very high accuracy results for one-dimensional convective simulation on practical (i.e., coarse) 
grids. The method automatically produces tight nonoscillatory resolution of discontinuities 
without distorting smooth profiles or clipping very narrow extrema. All these desirable fea- 
tures are obtained for little more cost than that of the basic third-order upwind scheme on 
which the algorithm is based (reference 2). Figure 1 shows what can be achieved in one dimen- 
sion. The figure shows computed points together with the exact solution for four profiles: a unit 
step, an isolated sine-squared wave 20Ax wide, a semi-ellipse 20Ax wide, and a narrow Gaus- 
sian (a = 1.94 Ax); at the time shown, the profiles have been translated by pure convection 45 
mesh-widths to the right (100 time-steps at a Courant number of 0.45). This particular example 
is based on third-order upwinding in smooth regions with automatic adaptive stencil expansion 
to seventh- or ninth-order upwinding locally, as needed, determined by monitoring the local 
absolute gradient across face ( i + ^), say, 

GRAD = |4> ;+1 - 4>j I (1) 

and the corresponding local absolute average curvature 

CURVAV = 0.5|(4, ( + 2 -<J> ( + ] )-(<l) [ -4>._ 1 )| (2) 

The algorithm also includes an automatic discriminator which decides when to apply the 
universal limiter (described in reference 3) and when to relax the limiter constraints. Ideally, 
the discriminator should activate the limiter in order to suppress unphysical numerical oscil- 
lations which would otherwise occur near sudden changes in gradient - each side of the step and 
at the "feet” of the semi-ellipse - without concomitant loss of resolution of the physical extrema 
(especially the narrow peak of the Gaussian profile). Clearly from Figure 1, this has been 
achieved; simulation of the sine-squared and Gaussian profiles is essentially exact, whereas 
resolution of the large-gradient regions of the other profiles is very tight. 
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Figure 1. One-dimensional pure convection of four profiles using the cost-effective 
3/7/9th-order scheme described in reference 1. 


The main purpose of the present paper is to show how to extend these desirable features to 
multidimensional flow problems. As is well known, one-dimensional algorithms do not auto- 
matically generalize to two and three dimensions simply by using the one-dimensional scheme 
component- wise in each direction. However, by using the concept of vector transient interpola- 
tion modelling 

4>(x , / + A/ ) = (J)(x — vA/,0 (3) 

fundamentally multidimensional convection schemes can be generated with the same prop- 
erties as their one-dimensional counterparts. Very high accuracy (in both space and time) can 
be obtained in a simple single-time-step explicit update formulation by converting time evolu- 
tion into a spatial interpolation problem at the earlier time-level, as represented by Equation 
(3). The multidimensional algorithms will be demonstrated in two dimensions. Once this is 
done, it becomes clear how to extend to three dimensions. To paraphrase a well-known aero- 
space quotation: algorithmically, there is a giant leap between one-dimension and two, but 
only a small step between two and three. 

In order to demonstrate the process of converting Equation (3) into conservative control- 
volume form, two-dimensional transient interpolation modelling will be considered for first- 
order upwinding, and three second-order schemes: Lax-Wendroff (reference 4), second-order 
upwinding and Fromm’s method (reference 5). The two-dimensional extension of the QUICK- 
EST scheme (reference 2) represents a uniformly third-order polynomial interpolation algo- 
rithm (UTOPIA). As with QUICKEST in one dimension, UTOPIA is susceptible to unphysical 
overshoots and undershoots if sudden changes in gradient are involved. Essentially nonoscil- 
latory results can be obtained by applying the universal limiter of reference 3 to the individual 
control-volume fluxes. Although UTOPIA has excellent phase accuracy, short-wavelength 
resolution is, of course, limited to third order. In principle, arbitrarily higher order resolution 
can be obtained locally - as in the one-dimensional scheme demonstrated in Figure 1 - by local 
adaptive stencil expansion. It appears that stencil expansion in a direction normal to a par- 
ticular control-volume face is much more effective than expansion in the transverse direction. 
Thus, although transverse terms are included to third order, higher order stencil expansion is 
taken only in the normal direction. Results for a third/seventh-order scheme are given for the 
well-known rotating-velocity-field benchmark test problem using three test profiles: a cylin- 
der, a cone, and a narrow Gaussian. In the results given here, an ad hoc discriminator is used in 
order to relax the limiter constraints near physical extrema. An automatic multidimensional 
discriminator is currently under development. 


228 




CONTROL- VOLUME FORM ULATION 

For simplicity, consider a two-dimensional square mesh (of unit grid size) with a control- 
volume using standard compass-point labelling. Equation (3) can be rewritten as 

<j)" + 1 = 4,(0,0,< + At) = - vAt, t) = $ n (-uAt, - vAt) ( 4 ) 

where <J> n (x,y) represents the behaviour of the convected field variable in the (upstream-biassed) 
vicinity of the control volume. The following algorithms depend on the functional form asumed 
for fy n (x,y). Assume that u and v are both positive and (for the moment) constant. 

First-Order Upwinding 

For example, consider the bilinear expression 


(}) n (x,y) = a + bx + ey + dxy (5) 

The four constants, a, b, c, and d, need to be evaluated; it is appropriate to use collocation at four 
nodal points with an upstream bias - in this case ( u and v positive), these would be 


K’ $w’ 4>«, and $ 


sw 


(6) 


as shown in Figure 2, giving 


4> n = 4 >£ + (4 >p -4>^)* + (4>p -4>g)y + 1(4 -4>^) - (4>$ - 4>s W ) l*y (7) 

as can be easily checked by putting x and y equal to 0 or -1, independently. Using Equation (4), 
the explicit update algorithm becomes (for u> 0 and v > 0) 


4> 


tl + 1 


4> 


c Xp 


-4V) 


C ,K 


-4>q) + 


; 4 [( 4 , p-4 , ^ ) - ( 4>s- 4>; w ) l 


(8) 


where c and c are the respective component Courant numbers. Equation (8) can be written in 
the general conservative control-volume form 


4> 


n + 1 


4> 


D + c 4> 

P X^W 


cd> -f c 6 — c <b 
xe y 1 s y n 


where the lower-case subscripts refer to face values, and 


(9) 


4 > e (j,j) = 4> U) (i+ Lj) 

and 


( 10 ) 


4>„ (». J) = 4> s (iJ+U 

guaranteeing convective conservation. This is achieved by identifying the face values as 


(ID 


4» 




(12) 
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Figure 2. Compass-point notation, showing nodes involved in the first- 
order upwind scheme (solid dots) for u,v > 0. 


and 


** = ~ f ( *s " ^ (13) 

again for u, v >0, of course. Equation (12), for example, can be rewritten in a form valid for 
positive and negative velocities as 


SGN(c ) c 

<b = <b L - 1 GRADN - - GRADT 

w w 2 2 

introducing the linearly interpolated face-value 

= \ «£ + 

the normal gradient across the west face 

GRADN = <t>" - 4>^ 

with upwind bias determined by the sign of c x 

SGN(c) = ± 1 fore >0 

X X ^ 

and the upwind-biassed transverse gradient 

GRADT = <b.., — d> m „ for c >0 and c >0 
Vv SW x y 

= - (}),„ for c >0 and c <0 

T NW T vv x y 

= d>„ — d>^ for c <0 and c >0 
x y 

= 4>„r “ for C <0 an d C <0 

N r x y 


(14) 


(15) 

(16) 

( 17 ) 

(18) 

(19) 

( 20 ) 
( 21 ) 
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A similar expression can be written for $ s . Then the east and north face values are obtained 
from Equations (10) and (11). Finally, if c and c appearing in Equation (14) (and the analo- 
gous expression for <J> S ) are interpreted as local face values of the respective component Courant 
numbers, the formulas can be considered to be valid for a spatially varying convecting velocity 
field. 

Second-Order Central 
Consider the second-order expression 

= a + bx + cx 2 + dy + ey 2 + fxy (22) 

The six constants are determined by collocation at six nodal points: and the four surround- 

ing points ($ N , $ s , and $ w ) together with one additional point. As with first order, the latter 
point is chosen on the basis of upwind bias: <Jwfor u>0, v>0; <Jw,for u> 0, u<0; <j> s/s for u< 0, 
u>0; and for u<0, u<0. The stencil is sketched in Figure 3(a) for u,v> 0. After evaluating 
the constants, using Equation (4), and rewriting in conservative control-volume form, the fol- 
lowing formula results for the west face 


c c 

<b = <b L - - GRADN - - GRADT 

w w 2 2 ) 

with a similar expression for <|> . As in the one-dimensional case, the difference between this 
and first-order upwinding, Equation (14), is the appearance of c x itself rather than SGNlc^) in 
the coefficient of the normal gradient term. Note that the transverse gradient term retains the 
same form. 

Equation (23) is the basis of the single-time-step explicit form of the Lax-Wendroff 
(reference 4) or Leith (reference 6) scheme extended to two dimensions in conservative control- 
volume form. The first two terms on the right of Equation (23) represent the one-dimensional 
formula; the transverse gradient term is the significant addition for two dimensions. It should 
be clear that in three dimensions there would be an additional upwind-biassed transverse grad- 
ient term in the z-direction (multiplied by —cJ2). The same would apply in the case of first- 
order upwinding. 

Other Second-Order Schemes 

If the stencil shown in Figure 3(b) is used to evaluate the constants in Equation (22), the two- 
dimensional form of second-order upwinding results. The one-dimensional form was originally 
discussed (in passing) by Fromm (reference 5) and was made popular in the aerospace industry 
by Warming and Beam (reference 7). The resulting formula (for the west face, for example) can 
be written 


c (1 — c ) c 

<t> = - - GRADN - CURVN - - GRADT 

<" «" 2 2 2 

where CURVN is the upwind-biassed normal curvature given by 


CURVN = <}>" - 24>; + for c * > 0 

or 

CURVN = 4>" - 24>" + 4>'^ for c^<0 


(24) 


(25) 

(26) 
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Figure 3. Second-order two-dimensional stencils (u,u>0). (a) Lax-Wendroff. 

(b) Second-order upwinding. (c) Fromm’s method. 

at the west face. Clearly, the appearance of the normal curvature term represents the 
significant difference between second-order upwinding and the second-order central formula. 
The same is true in one dimension. Fromm’s so-called "zero-average-phase-error” method 
(reference 5) represents a simple average between the second-order central and second-order 
upwind schemes. By comparing Equations (23) and (24), it is seen that this merely reduces the 
CURVN coefficient by an additional factor of 2. The stencil is shown in Figure 3(c). 

Third-Order Upwinding - UTOPIA 

In one dimension, Fromm’s method was an attempt to offset the lagging dispersion (phase error) 
of the Lax-Wendroff scheme by averaging it with second-order upwinding (containing inherent 
leading phase error). This was only partially successful; however, using the same one-dimen- 
sional stencil, it is possible to eliminate entirely the troublesome (third-derivative) dispersion 
term in the truncation error of the second-order schemes. The resulting explicit third-order 
upwind (QUICKEST) scheme has excellent phase behaviour; leading phase error stems from a 
small fifth-derivative term - and this is inherently damped by a fourth-derivative dissipation 
term (without introducing an artificial second-derivative diffusion term). The corresponding 
two-dimensional scheme is based on the third-order polynomial expression 

= a + bx + cx 2 + dx 3 + ey + f y 2 + g;y 3 + hxy + ix 2 ;y + jxy 2 (27) 

requiring 10 (upwind-biassed) collocation points. The appropriate stencil (for u,v >0) is 
sketched in Figure 4(a). The resulting formula for the west face value can then be written 


c (1 — c 2 ) c 

4> = <1> L - - GRADN - h CURVN - - GRADT 

w w 2 6 2 

2 

/ C C \ / C C C \ 

+ (---) CURVT + — - - TWIST 
V 6 4 / V 3 4 J 


where two new terms are evident. The upwind-biassed transverse curvature is given (for the 
west face) by 

CURVT = - 24>; + ^ ' for > 0 (29) 


CURVT = 4>" - 2(J>* + <t>" for c <0 

T /V r Y S x 
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Figure 4. UTOPIA stencils, (a) Complete stencil for u,u>0. (b) West face stencil for 
arbitrary u,v. (c) South face stencil for arbitrary u,v. 


The upwind-biassed "twist” term is given by 

TWIST = (4>" - 4$ - (<*>; - $ n sw ) for c > 0 (31) 

or 

TWIST = (<j£ - 4>p - (4>J W - 4>;) for < 0 (32) 

Figure 4(b) shows the stencil involved in computing the west face value when both positive and 
negative convecting velocity components are allowed for; Figure 4(c) shows the corresponding 
south face stencil. Extension to three dimensions requires additional GRADT, CURVT, and 
TWIST terms in an obvious manner. 


UNIVERSAL LIMITER 

As in the one-dimensional case, use of the two-dimensional second- and third-order 
schemes may result in unphysical oscillatory solutions if sharp changes in gradient are 
involved. In the case of the third-order scheme, this usually involves only relatively small (up 
to about 5%) overshoots or undershoots, at worst (reference 2). Even so, it seems desirable to 
eliminate this type of error (because of the possibility of nonlinear feedback instabilities in 
coupled equations). A universal limiting procedure, described in reference 3 for the one- 
dimensional case can be directly applied at each control-volume face. One first computes the 
multidimensional-stencil convected face value, given by Equation (28), for example. Then the 
"normalized” value is computed; for c x > 0 at the west face, this would be 





^LV 

K - O 


At the same time, the normalized adjacent upstream node value is computed 


(33) 


(<iV - Kw ] 
K - *ww' 


Then, if 0 ^ = 1, the normalized face value is constrained by 


(34) 





( 35 ) 
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and 


^ . - 
<j> < min(4> / c , 1) 

w W ‘ x' 


( 36 ) 


Outside of the monotonic range (i.e., if <|> w < 0 or > 1) one could use any simple nonoscillatory 
scheme such as 


rv rv> 



(37) 


Then the (unnormalized) face value is reconstructed using 


T£c x < 0, the normalized west face value is given by 

(4> — 4»P 

tv/ T UJ CL 

4 > = 

- V 

and the normalized adjacent upstream value is then 


(38) 


(39) 


4> t 


K 




4> 


E' 


(40) 


A similar procedure is used for the south face (based on the sign of c at that face). The east and 
north (limited) values are then given by conservation, as usual. 


ADAPTIVE STENCIL EXPANSION 

The uniformly third order polynomial interpolation algorithm (UTOPIA) described above 
is clearly limited in terms of short-wavelength resolution. In order to gain higher resolution, 
the same cost-effective strategy of adaptive stencil expansion - as used with such success in the 
one-dimensional case (reference 1) - can be used in two and three dimensions, as well. Stencil 
expansion in a direction normal to a particular control- volume face is fundamental for higher 
order resolution. Higher order transverse, twist, and other cross-terms, beyond the third-order 
terms of Equation (28), appear to have very little effect. Accordingly, suggested stencils for 
two-dimensional fifth- and seventh-order upwinding are shown in Figures 5(a) and 5(b), respec- 
tively (for u,v> 0). In other words, these are higher order one-dimensional schemes (reference 
3) applied component-wise, together with the complete third-order cross-terms; omission of the 
latter terms causes severe anisotropic distortion and significant loss of accuracy in velocity 
fields oblique (or skew) to the grid. 

Because of the (component-wise) one-dimensional nature of the proposed stencil expan- 
sion, the process can be automated by exactly the same procedure as used in the one-dimen- 
sional case. This is described in detail in reference 1. In the multidimensional code the one- 
dimensional adaptive stencil expansion criteria are applied independently at each of the west, 
south (and, in 3D, bottom) faces. In "smooth” regions, the respective values of GRAD and 
CURVAV are well below (pre-assigned) thresholds so that the (unlimited) UTOPIA scheme is 
being used. For most flows of practical interest, this will account for the overwhelming bulk of 
grid points, especially in three dimensions. Near isolated regions involving large values of 
GRAD or CURVAV at particular control-volume faces, thresholds will be exceeded, automati- 
cally switching the algorithm to fifth or seventh (or, in principle, arbitrarily higher) order at 
those particular points. 
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Figure 5. Stencil expansion (for u,v> 0) from fully third-order upwind 
to (a) quasi-fifth-order upwind; (b) quasi-seventh-order upwind. 

Two things should be noted. First, the higher order stencils will be needed only in very 
narrow regions (by definition of sharp change in value or gradient), thus requiring the more 
expensive computation at a relatively small number of grid points. This is an extremely cost- 
effective way to achieve very high accuracy on practical (i.e., coarse) grids - obviously an 
important consideration for three-dimensional simulation. Secondly, it should be clear that the 
location of the higher order stencils automatically changes as the flow evolves. As in the one- 
dimensional simulations described in reference 1, phase-accuracy is therey extremely tight. 
This is a critical attribute for any code designed to be used in large-eddy or full Navier-Stokes 
simulations. 


G-EXPANSION TECHNIQUE 

It is an instructive exercise to make a Fourier-von Neumann analysis (reference 8) of the 
multidimensional formulas discussed above and to compare Taylor expansions of their complex 
amplitude ratio (sometimes called "amplification factor”), G, with that of the exact solution (for 
constant v). This also gives some indication of how to incorporate higher order diffusion terms. 
For example, consider the exact complex amplitude ratio (reference 9) for the two-dimensional 
constant-coefficient convection-diffusion equation 

G= exp [— a (0 2 + 0 2 )] exp [ — i(c 0 + c 0 )] (41) 

x y x x y y 

where a is the nondimensional diffusion parameter and the 0’s are nondimensional wave- 
number components. A Taylor expansion of Equation (41) gives 

(c 0 + c 0 ) 2 (c 0 + c 0 ) 4 

G = 1 — a (0 2 + 0 2 ) — * + -- * 

* y 2 24 

(c 0 + c 0 ) 2 a (0 2 + 0 2 ) a 2 (0 2 + 0 2 ) 2 

x x y y x y x y 

+ 9 + 1) + — (continued)... 
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I (c 0 -f c 0 ) 

xx y y 


(c e + c e ) 2 ( C 0 + c e ) 4 

1 _ a(0 2 + 0 2 ) _ + 

x y 6 120 


(c 0 +c0 ) 2 a(0 2 + 0 2 ) a 2 (0 2 + 0 2 ) 2 

x x y y x y x y 

+ + — + 


( 42 ) 


At second order, note the appearance of the cross-term, c x Q x c 0 . This is accounted for by the 
GRADT term in the two-dimensional second-order formulas. y Component-wise application of 
one-dimensional second-order algorithms would miss this important term. Similarly, the cross- 
terms in the UTOPIA scheme are responsible for matching all convection terms in Equation 
(42) through third order in 0 x and 0 . It is a relatively simple matter to match diffusion terms to 
this order, as well, building bn experience with the one-dimensional QUICKEST scheme devel- 
oped in reference 2 and more fully explored in reference 10. In this case, the west face value, for 
example becomes 


c 

<j> = cJ) L - - GRADN 

W W 2 


d-cf) 


CURVN - - GRADT 
2 


+ (— — — + 


a \ / c c 

- ) CURVT + ( — 
2 / \ 3 


jJ TWIST 


(43) 


which should be compared with the pure-convection formula, Equation (28). The corresponding 
west face gradient is given by 


A . ^ 

Ax( — 
dx / w 


C C 

GRADN - - CURVN - - TWIST 


2 2 

The three-dimensional extensions of Equations (43) and (44) should, by now, be clear. 


(44) 


BENCHMARK TEST PROBLEMS 

The well-known "rotating-velocity-field” convection problem is used as a benchmark test. 
The velocity field is that of solid-body rotation so that a given initial profile should be swept 
around as if it were imbedded in a rotating solid. For pure convection, the exact solution is thus 
known. The following three initial profiles are considered: a cylinder with a base diameter of 
16 mesh-widths; a cone with the same base diameter; and a relatively narrow Gaussian distri- 
bution (o = 2 mesh- widths). The computation is carried out on a 55 X 55 grid with a maximum 
Courant number near 0.8. Figure 6 shows results of the Lax-Wendroff simulation after one-half 
rotation in the counter-clockwise direction; the exact solutions have been juxtaposed, for refer- 
ence. Note the typical trailing (phase-lag) oscillations, especially in the case of the cylinder. 
Figure 7 shows the corresponding two-dimensional second-order upwind simulation. In this 
case, phase-lead dispersion (partially obscured) occurs ahead of the simulated profile. Two 
nonoscillatory (TVD) schemes designed by Roe (reference 11) are shown in Figures 8 and 9. As 
seen, the Minmod results are quite diffusive for all profiles. Superbee does a reasonably good 
job on the cylinder, but tends to steepen and clip the other profiles. This is a well-known short- 
coming of second-order-based TVD schemes, especially those of supercompressive type that rely 
on negative artificial diffusion to enhance discontinuity resolution; this is explained in detail in 
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reference 3, Finally, figure 10 shows results that can be obtained using methods described in 
this paper. The results shown are for a seventh-order upwind scheme including all third-order 
(but no higher) cross-terms. An ad hoc discriminator is used to relax limiter constraints in the 
vicinity of physical maxima; as mentioned before, an automatic multidimensional discrimi- 
nator (similar to that described in reference 1 for one dimensional flow) is under development. 
This, of course, will be necessary before the code is applicable to general flow problems. Note 
from Figure 10 that resolution of the cylinder is better than that of Superbee, but without gross 
distortion of the other profiles. 


CONCLUSION 

The high-convection code described in this paper includes a number of features that are 
important for cost-effective accurate simulation of multidimensional unsteady flows on practi- 
cal grids. Being based on third-order upwinding, the code is totally free of artificial numerical 
diffusion (or viscosity); this is important because schemes based on artificial-viscosity methods 
are often solving the wrong problem - i.e., an artificially low- convection problem rather than 
the physical high-convection problem. UTOPIA contains all necessary cross-terms to third 
order, thus matching all terms in the Taylor expansion of the complex amplitude ratio through 
to third oi'der in the wave-number components. This guarantees excellent phase behaviour and 
isotropy regardless of the stream-to-grid angle. Codes based on second-order (or even fourth- 
order) central schemes inherently contain serious dispersion errors. Inclusion of the GRADT 
term should improve isotropy (in theory) - but this is usually masked by gross dispersion. 

Although far more accurate than first- and second-order methods, third-order upwinding 
may give rise to slight overshoots or undershoots. This tendency can be eliminated by using the 
universal limiter, developed in reference 3, on each control-volume face, independently. If 
higher order resolution is required, the strategy of adaptive stencil expansion - increasing accu- 
racy (above third order) only where needed - is extremely cost-effective. This is controlled by 
monitoring the absolute normal gradient and curvature across control-volume faces; as certain 
(pre-assigned) thresholds are exceeded, the code automatically switches to (in principle, arbi- 
trarily) higher order accuracy at the face in question. Finally, in order to give full resolution to 
local extrema, the limiter constraints need to be automatically relaxed in such regions (in addi- 
tion to using a higher order stencil). A fully automatic pattern-recognition discriminator of this 
type has been designed for one dimension (reference 1); the same principles appear to be appli- 
cable to multidimensional flows, as well, but a completely automatic multidimensional discrim- 
inator of this type is still under development. 
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Figure 6. Lax-Wendroff simulation after one-half rotation 
counterclockwise compared with exact result. 
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